Renormalization Group Approach to Non— equilibrium Green Functions in Correlated 

Impurity Systems 



T. A. Cost! 

Institut Laue-Langevin, B.P.156 38042 Grenoble, Cedex 9, France. 
Universitat Karlsruhe, Institut fur Theorie der Kondensierten Materie, 76128 Karlsruhe, Germany. 



On 
ON 

3 



i 

c3 



o 
o 



> 

on 
o 

NO 

o 

On 



c 

o 
o 



X 
S3 



We present a new technique for calculating non- 
equilibrium Green functions for impurity systems with local 
interactions. We use an analogy to the calculation of response 
functions in the x-ray problem. The initial state and the final 
state problems, which correspond to the situations before and 
after the disturbance (an electric or magnetic field, for exam- 
ple) is suddenly switched on, are solved with the aid of Wil- 
son's momentum shell renormalization group. The method 
is illustrated by calculating the non-equilibrium dynamics of 
the ohmic two-state problem. 

PACS numbers: 71.27.+a,73.20.Dx,71.10.+x,72.15.Qm 



I. INTRODUCTION 

Recently there has been interest in the non- 
equilibrium transport properties of small devices, such 
as quantum dots and ultrasmall tunnel junctions KB- 
These systems, together with some resonant tunneling 
devices 0, offer new possibilities for studying many- 
body effects due to strong local Coulomb interactions. 
The importance of these interactions in small devices 
is seen, for example, in the suppression of tunneling or 
Coulomb Blockade effect p|. Phenomena such as the 
Kondo effect in quantum dots and the Fermi edge singu- 
larity in resonant tunneling devices have been predicted 
P-|TT| and some aspects of these have been confirmed 
experimentally 0. The usual starting point for dealing 
with non- equilibrium transport in such systems has been 
the formalism developed by Keldysh jl2) and Kadanoff 
and Baym [jl3|. Below we describe a non-perturbative 
approach based on the numerical renormalization group 
(NRG) method JlJ,|l^] which allows the calculation of 
non-equilibrium Green functions for the above systems. 
We consider only the case in which the perturbation (an 
electric or magnetic field) causing the non-equilibrium 
effects is suddenly switched on at time t = 0. The non- 
equilibrium Green Functions will then be calculated by 
solving an initial (t < 0) and final (t > 0) state prob- 
lem as in the case of calculating the photoemission and 
absorption spectra in the x-ray problem jl6| . In this pa- 
per we concentrate on calculating the non-equilibrium 
properties of a specific model, the Ohmic two state sys- 
tem JIt]]. The application of the method to the systems 
mentioned above follows along the same lines, the only 
difference being the solution, using the NRG method, of 



different initial and final state Hamiltonians. 

The paper is organized as follows: in Sec. II we intro- 
duce the standard model of the Ohmic two state system, 
formulate the problem of calculating the non-equilibrium 
dynamics of this model in terms of solving an initial 
and a final state problem and introduce an equivalent 
model, the anisotropic Kondo model, which we actually 
use in the calculations. In Sec. Ill we describe the NRG, 
its application to dynamical quantities and an approxi- 
mate evaluation of the formally exact expressions for the 
non-equilibrium quantities. An exact evaluation of non- 
equilibrium quantities first has to overcome certain tech- 
nical difficulties, which we describe, and is postponed 
for the future. Sec. IV contains our results for the 
non-equilibrium dynamics of the Ohmic two-state sys- 
tem, obtained on the basis of NRG calculations for the 
anisotropic Kondo model. In Sec. V we summarize our 
main conclusions. 



II. FORMULATION 

A. The Ohmic Two-State System 

The non-equilibrium properties of the Ohmic two state 
problem are of main interest in macroscopic quantum 
coherence experiments in rf SQUIDs |fl8|| . Typically, an 
rf SQUID can be in one of two possible ffuxoid states. By 
applying a bias (corresponding to an external magnetic 
field), for times t < 0, the system is prepared in one of 
the two states. The dynamics after the bias is removed 
at t > is then intrinsically a non-equilibrium property. 
The Hamiltonian, H , of the system is time dependent 
with a sudden perturbation at t — 0, so that we can 
write H{t) = [1 - 9(t)]Hj + 6(t)H F where Ht,H f are 
the Hamiltonians before and after the bias is switched 
off. The Hamiltonian Hj p describing the Ohmic two 
state system is given by the spin-boson model JTsf ] , 
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Here <ji,i — x,y, z are Pauli spin matrices, the two states 
of the system correspond to a z =\ and o z =\ (i.e. 
a z =1,1 correspond to the two possible ffuxoid states 
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of the rf SQUID). A is the bare tunneling matrix ele- 
ment and e is a bias. The environment is represented 
by an infinite set of harmonic oscillators (labeled by the 
index a) with masses m a and frequency spectrum u a 
coupling linearly to the coordinate Q = \q$(J z of the 
two-level system via a term characterized by the cou- 
plings C a (the two-level system coordinate could be the 
total flux, <fi — <fii(T z , in the case of an rf SQUID ex- 
periment). The environment spectral function is given 

in terms of these couplings, oscillator masses and fre- 

c 2 

quencies by J{uj) = § I] Q (^Sr)^( w ~ ^a)- In the case 
of an Ohmic heat bath, of interest to us here, we have 
J(lu) = 2naLU, for lu << lo c , where uj c is a high en- 
ergy cut-off and a is a dimensionless parameter char- 
acterizing the strength of the dissipation. This form for 
the spectral function is appropriate for describing quan- 
tum dissipation in an rf SQUID. Preparation of the sys- 
tem in a state with a z = +1 with the oscillators re- 
laxed about this state is equivalent to setting e = — oo in 
Hsb, so the initial state problem corresponds to solving 
Hi = Hsb{£ — — oo). Similarly, the final state problem 
is given by Hp = Hsb(c = 0). The Ohmic spin-boson 
model has been intensively studied (for reviews we refer 
the reader to [^9|j2C| l). We outline some of its features in 
order to introduce some useful notation. The model has a 
low energy scale, A r < A for A << u> c , which depends on 
the dissipation strength a, and which may be interpreted 
as a renormalized tunneling amplitude. For a << 1 
the dynamics corresponds to damped oscillations, with 
a crossover to incoherent behaviour with increasing dis- 
sipation strength. For a — > 1~, the renormalized tun- 
neling amplitude vanishes giving rise to the phenomenon 
of "localization" or "self-trapping" for a > a c ~ 1 (a c 
depends also on the precise value of A). The dynami- 
cal quantities exhibiting the above features are defined 
below. 

B. Non— equilibrium Dynamical Quantities 

The simplest non-equilibrium dynamical quantity to 
study for the spin-boson model is the quantity P(t) = 
{&z{t)) Pl Jl9| where the thermodynamic average is taken 
with respect to the initial density matrix pi(fl) = 
e -@Hi jTr e~P Hl , (3 is the inverse temperature, and the 
time evolution is with respect to the Hamiltonian af- 
ter the bias is switched off at time t = 0, i.e. cr z {t) — 
e lHFt a z e- tHFt . Hence, P(t) = 1 for t < due to the infi- 
nite bias e = — oo, and for t > 0, when the bias is switched 
off (e = 0), P{t) describes how the two-level system co- 
ordinate <7 2 relaxes to its long-time value of zero. An- 
other quantity of interest is the retarded two time Green 
function, G T (t,H) = -iO(t - f)([a z (t),a z (t')]) pi , with 
the thermodynamic average defined as above. Since time 
translational invariance is broken, G r (t,t') depends on 



both times explicitly. We consider the Fourier transform 
of G r (t, t') with respect to both the sum t + t' and differ- 
ence t — t'oi the time variables. The resulting spectral 
density C r (uj,D,) = — ^Im G r (ui + iS,il), with uj,n, the 
Fourier frequency variables corresponding to t—t', t+t', is 
given within a Lchmann representation by the following 
expression, 

E — E a 

x (■m F \a z \m' F )(m' F \a z \mp)S(n — — — ) 

x [S (u;+[ EmF+ 2 Em - -E m , F \) 

-6(.-[ EmF+ 2 E ^ -E m i p })}. (2) 

Here, E mi and \mi) are the many-body eigenvalues and 
eigenstates of the initial state Hamiltonian Hj, Zj the 
corresponding partition function and E mF , \mp), E rn ' F , 
|m' F ), . . . the many-body eigenvalues and eigenstates of 
the final state Hamiltonian Hp. In the equilibrium case, 
Hj = Hp = H, and the corresponding spectral density 
C^ q (iv,fl) reduces to, 

C e r «(u,Sl) = y E e^ E -\(m\a z \m')\ 2 S(n) 

m.m' 

x [6(u + [E m - E m i}) - 5{lo - [E m - E m i})}, (3) 

where E m , |m) are the many-body eigenvalues and eigen- 
functions of H and Z is the corresponding partition func- 
tion (the delta function, 5(Q), in the above expression 
reflects the fact that in the equilibrium case, G r (t, t') de- 
pends only on the difference of the time variables). 

We see that the non-equilibrium spectral density dif- 
fers from the equilibrium one in several ways: first, even 
at T = 0, no groundstate energy appears in the delta 
functions, the excitations are between arbitrary (final) 
excited states of the system. This reflects the fact that 
there is no stationary groundstate for a non-equilibrium 
situation. Secondly, the non-equilibrium aspects, which 
are a result of an initial state preparation, are reflected 
in the presence of overlap matrix elements between the 
initial and final states. Finally the dependence on f2 is 
a measure of the importance of transient effects. Ne- 
glecting these effects results in the following simplified 
expression for the spectral density, 

C 0r (Lu) = C r (u,Q = 0) 

= Y E e-^\(mi\mp)\ 2 \(mp\a z \m' F }\ 2 S(Q) 

x [5(cj + [E mF - E m i F ]) - 5(lu - [E mF - E m i F ])). (4) 

This describes the steady state case t + t' — > oo. In a 
strict sense it is not a non-equilibrium quantity, although 
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it does take into account the effects of an initial state 
preparation on the correlation function (a z (t)a z (0)) . Our 
motivation for calculating this quantity is simply to illus- 
trate that our technique applies also to two-time Green 
functions. The calculation of the full spectral density 
C r (u>,Q,), including both frequencies involves a straight- 
forward generalization. 

Similarly we can write the Fourier transform of P(t), 
within a Lehmann representation, as, 

1 



Y X^ c feT Cfe 'T ~ c k.i c k'l)S z + 9»BhS z 



(6) 



P Em ' (mi\m F )(m' F \mi) 
x (m F \a z \m' F )S(uj - (E mF - E m > )), (5) 



where the same notation as above is used. 

P(t) = J Q P(u>) cos(u)t)dt contains information on the 
onset of quantum oscillations in the two level system. 
For small values of the dissipation strength, a << 1, 
it is known from the "Non-Interacting Blip Approxima- 
tion" (NIB A) fH] that P(t) exhibits damped oscillations 
with a renormalizcd tunneling frequency A r = A[— ] 1 ~ a . 
P(u>) will exhibit two peaks at ui = ±A r . At the exactly 
solvable Toulouse point plj, a — 1/2, where the Ohmic 
two-state system reduces to the resonant level model 
p9|p2| , the dynamics is incoherent and P(t) decays ex- 
ponentially. P{u>) consists of a single peak at u> = 0. It 
is not clear at which value of the dissipation strength the 
crossover to incoherent behaviour occurs, in particular 
whether it occurs at exactly a = 1/2 or for some smaller 
value of a. This may depend on the definition of the 
crossover and on whether equilibrium or non-equilibrium 
quantities are being studied. For equilibrium quantities 
a smooth crossover has been found to occur at a = 1/3 



C. The Anisotropic Kondo Model 
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The first term represents non-interacting conduction 
electrons and the second and third terms represent an 
exchange interaction between a localized spin 1/2 and 
the conduction electrons with strength J±,Jn. A local 
magnetic field, h, coupling only to the impurity spin in 
the Kondo model (the last term in Eq. ^) corresponds to 
a finite bias, e, in the spin-boson model. The correspon- 
dence between H and H$b is then given by e = gpBh, 



S is 



= p,Jj_ and a = (1 + ^) , where tan S 
the phase shift for scattering of electrons from a potential 
Jy/4 and p — 1/2D is the conduction electron density of 
states per spin at the Fermi level for a flat band of width 
2D 19 We note that weak dissipation (a — > 0) in 



the spin-boson model corresponds to extreme anisotropy 
(J|| — > oo) in the Kondo model. For zero dissipation, 
J,, = oo, the two states V± = ^(1 T)l Do ± II) | T)o) 
made up from the impurity states and the local con- 
duction electron Wannier orbitals |er)o = J2k c ka\ vac )> 
where \vac) is the vacuum, are split by Jj_ = A (with 
the identification lo c = 2D) and are completely decou- 
pled from the rest of the conduction band, thus forming 
an isolated two-level system. The system exhibits coher- 
ent oscillations with P(t) = cos(At). As Jy is decreased 
from +oo, the two-levels become weakly coupled, with 
strength -j- tx a, for pj^ >> 1, to the remaining con- 
duction states and their splitting is renormalized down- 
wards. The low energy scale of the model is given by 
the Kondo temperature for the anisotropic Kondo model, 
Tk(J±,J\\) < J_i_, for Jl << D, which in the lan- 
guage of the dissipative two-state system corresponds to 
a renormalizcd tunneling amplitude A r J2~5{ |. An exten- 
sive discussion of the equivalence between the anisotropic 
Kondo model and the Ohmic two-state system is given 
elsewhere. 



Instead of solving directly the spin-boson model with 
the NRG method it is more convenient to solve an equiva- 
lent fermionic model which has the same dynamics. This 
is the Anisotropic Kondo Model ( AKM) . The equivalence 
has been shown at the Hamiltonian level via bosonization 
|p3| . This was believed to be valid in the region a > 1/2, 
which corresponds (see below for the precise statement 
of the equivalence) to the region in the parameter space 
of the AKM between weak-coupling and the Toulouse 



point. In fact, recent work |24| shows that the equiva- 
lence extends beyond the Toulouse point into the region 
describing weak dissipation < a < 1/2 (or large anti- 
ferromagnetic J\\ in the AKM, see also Q). The AKM 
is given by Q 

fe.tr kk' 



III. CALCULATION OF NON-EQUILIBRIUM 
DYNAMICS VIA THE NRG 

A. The NRG 

The dynamical quantities P(t) and G r (t,t') defined 
above for the two-level system translate into the corre- 
sponding quantities for the Kondo model (a z — > S z under 
the equivalence). We calculate these quantities by ap- 
plying Wilson's momentum shell renormalization group 
method generalized to the calculation of dynamical quan- 
tities (e.g. |24|,|2j|). Thus in addition to solving an initial 
state problem Hj = Hakm(^ = — oo) and a final state 
problem Hp — Hakm{^ — 0), the final state matrix ele- 
ments of the variable a z and the overlap matrix elements 
appearing in the above expressions for P(u>),C r (u>) are 
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also calculated. The diagonalization of Harm proceeds 
as follows (for full details see (i) the spectrum is lin- 
earized about the Fermi energy — > vpk, (ii) a logarith- 
mic mesh of k points k n = A~™ is introduced to achieve 
a separation of energy scales, and (iii) a unitary trans- 
formation of the Cfc CT is made such that /o CT = Efc Cka 
is the first operator in a new basis, f n<T ,n = 0,1,..., 
which tridiagonalizes H c — J2k^ e fcM c l/j c fcp m k-space, 
i-e- H c -> E,EZo^A- n/2 (fl +1 ,fn» + h.c), with 
£„ -> (1 + A- 1 )/2 for n » 1. The Hamiltonian (§) 
with the above discretized form of the kinetic energy is 
now diagonalized by the following iterative process: (a) 
one defines a sequence of finite size Hamiltonians Hn = 

T,n=0 Cn^ n/2 (fl +lfl fn^ + h.C.) + ^(/ \/o^- + 

/oi/ot^ + ) + J i (/o T /oT - flihl)S z for N > 0; (b) the 

JV — 1 

Hamiltonians Hn are rescaled by At such that the 

— N — 1 

energy spacing remains the same, i.e. Hn = A^~Hn- 
This defines a renormalization group transformation 
H N +i = A l l 2 H N + ^ ^(/Ir+i^ + &.c.) - ^g,jv +1 , 
with Eg,n+i chosen so that the ground state energy of 
Hn+i is zero. Using this recurrence relation, the se- 
quence of Hamiltonians Hn for N = 0, 1, ... is iteratively 
diagonalized. This gives the excitations and many body 
eigenstates at a corresponding set of energy scales ljn de- 

N — 1 

fined by cjat = A and allows a direct calculation of 
the dynamical quantities P(t), G r (t, t'), or more precisely 
their Fourier transforms using the Lehmann representa- 
tions (||,||). Our results were obtained for A = 2, keeping 
the 320 lowest states at each iteration. Truncating the 
spectrum in this way restricts the range of excitations lj 
at iteration N to be such that wjy < w < Kljn, where 
K = if (A) w 7 for A = 2. In this paper we discuss only 
the T = results. 

In diagonalizing the Hamiltonians Hn the following 
symmetries are used to reduce the size of the matrices, 
(i) conservation of z-component of total spin = S z + 

E«=o1(/It /t - /Ii /l ) > ( u ) conservation of total pseudo- 
spin [|30f , where the pseudo-spin operators j + ,j~,j z are 
defined by j+ = E^oM)"/,^, j~ = U + V and 
3z = J2n=o( fmfniJ,- D- These symmetries hold for both 
the zero and finite bias cases, which we consider in this 
paper (e = for the final state problem and e = — oo 
for the initial state problem). At this point we men- 
tion that although it is in principal possible to apply this 
renormalization group method directly to the spin-boson 
model, in practice there are disadvantages in the case of 
the bosonic problem [3l] . 

B. Broadening procedure for the discrete spectra 

We note that the use of a discretized model implies 
that at each iteration P(u>) is a series of delta functions. 
Smooth curves are obtained by broadening these delta 



functions with a Gaussian of width appropriate to the 
level spacing of Hn p9[| . The width of the Gaussians 
will not influence the intrinsic peak widths at low ener- 
gies, where the logarithmic spacing ensures high resolu- 
tion, but may do so at higher energies where the resolu- 
tion is lower. This problem has been discussed in detail 
in p3[ , where a refined NRG scheme to overcome it has 
been suggested. In most problems one is interested in 
the low energy behaviour and such a refinement is not 
required, however for the spin-boson model considered 
here there is interesting dynamics at high energies and 
such a refinement is required in order to obtain a com- 
plete description of all aspects of the high energy (short 
time) dynamics (by high energy, we mean high relative 
to the renormalized tunneling amplitude, but still low 
relative to the high energy cut-off in the model) . In the 
present paper we have not implemented the above re- 
fined scheme, so we shall state below which aspects of 
the high energy dynamics we are able to obtain within 
the unmodified NRG scheme. 

C. Approximate evaluation of non— equilibrium 
quantities 

At this stage we point out an importance difference 
in the calculation of non-equilibrium quantities with re- 
spect to equilibrium quantities such as C^. q {oj, Q). At 
T = 0, the latter can be calculated on each energy scale 
lon by restricting attention to a single energy shell N. 
This is due to the existence of a stationary groundstate, 
from which all excitations in the expression for C^ q can 
be measured. The delta functions in the Lehmann rep- 
resentation, (^|), then imply that in order to calculate 
C^ q at frequency to, only one energy shell, that for which 
uj ss ujn, is required, i.e. 

C?MT = Q) = I^KGSIa.KM 2 ^) 

m' 

x [8(u> + [Eg s - - 5{u: - [E» s - E»])], (7) 

where |m')jV) \GS >n arc the eigenstates and the 
groundstate, corresponding to iteration N in the NRG 
procedure, E%, and Eq S the corresponding eigenvalues, 
and Zn the (T = 0) partition function. Contributions 
from energy shells n — 0,1, ... , N—l have E^ — E^g > lj 
so they make no contribution to the spectral density C^ q 
at frequency lj s» ljn- At finite temperatures, there can 
be contributions from higher energy states, although the 
Boltzmann factors in (^|) will suppress these. 

The situation is different for non-equilibrium dynam- 
ical quantities, such as P{ui). From the Lehmann rep- 
resentation (||) we see that, even at T = 0, no ground- 
state energy appears in the delta functions. Instead the 
(T = 0) excitations are between arbitrary (final) states. 
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The response at frequency lu ~ lon can have contribu- 
tions from all energy shells n = 0, 1, . . . , N. Hence, if we 

N — 1 

evaluate P(lo), at the frequency ui ~ ujjy = DA ~ , by 
taking into account only the iV'th energy shell, we have 
only the approximate result 



1 



JN,I 



E 



(™/,Gsl^F)Ar(m^|TO/ !GS )Ar 



x (m F \a z \m' F ) N 5{u-{E? nF -E")). 



(8) 



The approximation implicit in this procedure is that 
small excitations u> ~ i?JJ, F — £7™, between higher en- 
ergy excited states (i.e. from energy shells with n < N), 
make a negligible contribution to P(u>) compared to those 
on a scale u> ~ lom from the iV'th energy shell. This 
is clearly only an approximation to the formally exact 
expression (|^). The results and arguments to be pre- 
sented below show that the above approximation has a 
regime of validity and a regime where it breaks down. 
We stress, however, that a full multiple-shell evaluation 
of (||) will ultimately be required. This is currently not 
feasible within the standard NRG procedure described 
in Sec. Ill A-B. One first has to overcome two problems: 
(i) a possible double counting of excitations in adding 
contributions to P(lo) from different energy shells, (ii) 
a meaningful way of adding contributions to P(uj) from 
higher energy shells which have different resolutions. The 
second problem can be overcome by replacing the stan- 
dard NRG procedure of Sec. Ill A-B by one which elim- 
inates any dependence of static and dynamic quantities 
on the logarithmic discretization . These problems do 
not arise in the case of equilibrium dynamical quantities 
since, as discussed above, the spectral densities can be 
calculated essentially without approximation by restrict- 
ing attention to just one energy shell for each frequency of 
interest pjl (although, as discussed in Sec. Ill B, the use 
of a logarithmic discretization and a Gaussian broaden- 
ing procedure can overestimate the widths of high energy 
peaks in spectral densities). 

Returning to (||) we expect this to be a valid approx- 
imation as long as orthogonality effects do not signifi- 
cantly affect the overlap matrix elements (m'p\mi 1 Gs)N 
appearing in the above expression (and similar expres- 
sions for Co r (w)). When this occurs, as it will do for 
sufficiently low energies [p2[ , an increasing number of fi- 
nal states \rriF >n will be nearly orthogonal to the ini- 
tial groundstate \mj t GS >n and it will be necessary to 
include higher energy states. Information on higher en- 
ergy states is contained in previous iterations within the 
NRG procedure. Thus if we evaluate the non-equilibrium 
quantities approximately by using only one energy shell, 
then we will obtain results with a restricted range of va- 
lidity. The range of validity can be estimated: orthogo- 
nality effects between initial and final states become im- 



portant in the strong coupling regime, i.e. for frequen- 
cies u> « Tk, where Tk is the low energy scale of the 
AKM (or the renormalized tunneling amplitude, A r , in 
the language of the spin-boson model). However, the 
exponent with which the matrix elements (m' F \mi t Gs) n 
vanish will also influence the range of validity. We find, 
by keeping only one energy shell for each frequency u), 
that P(uj) ~ Cor(w)/w ~ \ui\ « A r (noticeable in 
the low energy part of our results in Fig. 1 and Fig. 2, 
presented in the next section). The vanishing of P{lo) 
with lu — > reflects the orthogonality of the matrix el- 
ements in the expression for P(lu) and we see that the 
exponent governing this is a/2. The exact behaviour of 
P(u>) and Co r (u>) for lu — > is not rigorously known. 
In special cases, such as at the Toulouse point, a = i, 



it is known that P(ui — > 0) ~ [ 



enst, and 



P(t) ~ e , t — + oo. This type of exponential relaxation 
is expected for other a in the range < a < 1. In any 
case, we see that the approximate evaluation of the non- 
equilibrium dynamics taking just one energy shell into 
account for each frequency is expected to be accurate for 
weak dissipation (when the overlap exponent is small) 
and for energies which are not too small relative to A r . 
In the context of macroscopic quantum coherence exper- 
iments in SQUIDs, one is interested in the short time 
dynamics for times up to approximately 1/A r . The long 
time behaviour t >> 1/A r is of interest in other contexts 
(e.g. microscopic two level systems), and for these it will 
be necessary to carry out the full calculation, including 
higher energy shells, for the non-equilibrium dynamics. 
We expect that such a calculation will give results as ac- 
curate as those for equilibrium quantities [E4| . 



IV. RESULTS 

For weak dissipation, a << I, one expects damped os- 
cillations of the two level system at a frequency reduced 
relative to the bare tunneling frequency, A, due to the 
coupling to the environment. From Fig. I, which shows 
P(u>), we see that this expected behaviour is reproduced 
by our method. The presence of an inelastic peak in 
P(uj) at io — A*, indicates damped oscillations of fre- 
quency A*. The width, 7*, of the inelastic peak gives 
the characteristic time I/7* for the decay of the oscil- 
lations. For a << I we find a renormalized tunneling 
frequency A* w A r < A with A r = AfA-jT?^ anc i A r 
is the relevant low energy scale for weak dissipation as 
obtained within the NIBA and several other approaches 
( p4| , |35|| ) . The larger renormalization of the tunneling 
amplitude with increasing dissipation is clearly seen in 
Fig. 1. Qualitatively similar results are found for the 
quantity Co r (uj)/io (shown in Fig. 2). For weak dissipa- 
tion, the inelastic peak width 7* ~ aA* vanishes linearly 
with a icil. 
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FIG. 1. P(u) for A = 0.1 (Jr. = 0.10397) a = 0.1 
(J|| = 4.698), solid line, and a = 0.2 (Jy = 3.12759), dashed 
line. Energies are in units of D = ^ = 1. The relation be- 
tween the AKM parameters J± \\ and the spin-boson model 
parameters a, A is given in the text. Renormalizations of 
these parameters due to the discretization of the AKM model 
have been described in p4| 

Consequently, for a << 1, the use of a logarithmic 
discretization does not provide the necessary resolution 
at the peak position which is needed for resolving the 
intrinsic peak width, i.e. any broadening of the discrete 
spectra will overestimate the inelastic peak width. This 
problem may be overcome by averaging over discretiza- 
tions using the procedure in f33|j . The standard proce- 
dure used here gives correctly the positions and weights 
of the delta functions in the discrete spectra (this has 
been shown in detail for all a in the range < a < 1 for 
the case of equilibrium dynamical quantities (24) ) . The 
decay of the weights of the delta functions with increas- 
ing energy and hence the frequency dependence of -P(w) 
at energies A r < oj < u> c is also correctly captured by our 



procedure. We find that P(u>) ~ Cbr(w)/w 



-(3-2a) 

for A << uj c , A r <C w <C lo c and for < a < 1 (the 
error in the exponent is typically less than 0.1%). Thus, 
the short time behaviour of P(t) is identical to the NIBA 
result P(t) = 1 - ct 2 ^-^ for u c r x < t < A r _1 . This 
gives independent confirmation that the NIBA is correct 
for P(t) at short times. 

Hence for weak dissipation we recover the known pic- 
ture Jl^,^(| of damped oscillations at reduced tunneling 
frequency A r . In particular, the short-time dynamics is 
non-trivial in the sense that P(t) and correlation func- 
tions depend on a dependent exponents. On increas- 
ing the dissipation strength, the inelastic peak in P(u>) 
narrows and the incoherent contribution (P(lu = 0)) be- 
comes larger. For sufficiently strong dissipation, we ex- 
pect the incoherent part of P(uj) to dominate and lead to 
incoherent dynamics of the two level system. Although 
the question of the crossover from coherent to incoher- 
ent dynamics in equilibrium quantities has been investi- 



gated with high accuracy using the NRG p4] , the same 
question is technically more difficult for non-equilibrium 
quantities such as P(oj). In this paper we have calcu- 
lated the non-equilibrium dynamics only approximately, 
and the approximation used, which we described in detail 
in Sec. Ill C, has a range of validity which restricts us to 
weak dissipation and energies which are comparable to or 
higher than the low energy scale of the model. To address 
the question of a crossover from coherent to incoherent 
dynamics we need to evaluate P(u>) taking into account 
several energy shells for each frequency lo, as detailed in 
Sec. Ill C. 
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FIG. 2. Cor(w) for A = 0.1 (J± = 0.10397). (a) a = 0.1 
(J,! = 4.698), solid line, and (b) a = 0.2 (Jy = 3.12759), 
dashed line. Energies are in units of D = ^ = 1. 



V. CONCLUSIONS 

To summarize, we have presented a numerical renor- 
malization group approach to the calculation of non- 
equilibrium Green functions in correlated impurity sys- 
tems. The approach uses an analogy to the calculation 
of response functions in the x-ray problem where the dis- 
turbance is sudden. The method was illustrated by calcu- 
lating the non-equilibrium dynamics of the Ohmic two- 
state system. An approximate evaluation of the non- 
equilibrium quantities, taking only one energy shell into 
account for each frequency range, gave accurate results 
within the expected range of validity of the approxima- 
tion: specifically the approximate evaluation of P(t) gave 
P{t) = 1 - ct 2 ^-^ for uj- 1 < t < A" 1 for all dissipa- 
tion strengths provided A << uj c . This is in agreement 
with the NIBA prediction, which is known to be accu- 
rate at short times. Including additional energy shells 
in the evaluation of non-equilibrium quantities should 
give essentially exact results. This is left for future work 
as it requires a modified NRG procedure (outlined in 
Sec. Ill C). The method is non-perturbative and can 



G 



be used to study the effects of local interactions on the 
non-equilibrium transport through small devices such as 
quantum dots and tunnel junctions. An important as- 
pect of the method is that the disturbance (an electric or 
magnetic field) can be arbitrarily large, so it may be used 
to study the non-linear regime in the I-V characteristics 
of small devices, such as quantum dots. 
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